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Abstract 

We  study  the  behavior  of  NURBS-based  Isogeometric  Analysis  on  problems  of  large-deformation  plasticity. 
We  evaluate  the  performance  of  standard  NURBS  elements  and  elements  based  on  the  F  formulation  of 
Elguedj  et  al.  (T.  Elguedj,  Y.  Bazilevs,  V.M.  Calo,  T.J.R.  Hughes,  B  and  F  projection  methods  for  nearly 
incompressible  linear  and  non-linear  elasticity  and  plasticity  based  on  higher-order  NURBS  elements,  Com¬ 
puter  Methods  in  Applied  Meehanies  and  Engineering^  197  (2008),  2732-2762).  We  determine  that  standard 
measures  of  evaluation  employed  in  the  literature,  namely,  displacements  at  selected  locations  and  graphs 
of  reaction  forces  versus  displacements,  are  often  misleading.  On  the  other  hand,  stress  distributions,  in 
the  form  of  contour  plots,  are  the  most  revealing  measure  of  element  performance.  We  also  determine  that 
the  concept  of  “mesh  locking” ,  which  has  dominated  investigations  of  low-order  elements,  is  not  a  relevant 
issue  for  higher-order  NURBS  elements  for  problems  of  large-deformation  plasticity.  However,  standard 
higher-order  NURBS  elements  of  type  Qk,  of  continuity  class  /c  >  2,  typically  exhibit  spurious  stress 

oscillations,  whereas  the  F  elements  of  type  Qk/Qk-1  produce  good  results  in  all  cases. 

Key  words:  Isogeometric  analysis,  NURBS,  multiplicative  plasticity,  projection  method,  F-bar 


1.  Introduction 

The  first  applications  of  IsoGeometric  Analy¬ 
sis  (IGA)  to  nearly-incompressible  elasticity  prob¬ 
lems  were  presented  in  Elguedj  et  al.  [1]  in  which 
displacement-based  B  and  F  formulations  for  small- 
and  large-deformation  problems  were  developed. 
These  methods  employed  a  lower-order  space  for 
representing  volumetric  deformations  in  order  to 
avoid  mesh  locking  due  to  nearly  isochoric  kine¬ 
matics.  The  B  formulation  developed  was  stan¬ 
dard,  following  the  original  ideas  of  Hughes  [2,  3]. 
However,  the  F  formulation  was  novel  and,  in  our 
opinion,  was  an  improvement  over  what  had  ex¬ 
isted  previously.  The  basis  of  its  derivation  was  a 
modified  minimum  potential  energy  principle,  the 
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modification  being  the  representation  of  the  defor¬ 
mation  gradient  in  terms  of  an  “improved”  F  .  The 
consistent  tangent  stiffness  matrix  was  shown  to 
be  symmetric  and  contained  several  “unexpected 
terms”  (see  Elguedj  et  al.  [1]  for  further  details). 
When  combined  with  a  Newton- Raphson  solution 
strategy,  it  exhibited  second-order  iterative  con¬ 
vergence,  as  theoretically  predicted.  The  selec¬ 
tion  of  the  spaces  of  basis  functions  was  based  on 
the  desire  to  include  and  generalize  the  most  uti¬ 
lized  finite  elements  for  nonlinear  solid  mechanics 
computations,  namely,  the  four- node  quadrilateral 
and  eight-node  hexahedral,  “constant  pressure  ele¬ 
ments”  ,  which  have  been  the  stalwarts  of  nonlinear 
commercial  and  research  codes  for  approximately 
40  years.  The  basis  functions  may  be  described  as 
“Qk/Qk-1”  that  is,  the  displacement  are  denoted 
Qk,  and  are  -continuous,  whereas  the  volu¬ 
metric  measures  of  strain,  namely,  the  dilatation 
in  the  small-deformation  case,  and  the  jacobian  of 
the  deformation  gradient  in  the  large-deformation 
case,  are  denoted  Qk-1  and  are  C^“^-continuous. 
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Clearly,  when  k  =  1  this  reduces  to  our  version  of 
the  constant  pressure  element,  that  is,  Ql/QO.  The 
higher-order  and  smooth  versions  of  these  elements, 
for  which  /c  >  2,  were  previously  unexplored,  but 
were  shown  to  perform  very  effectively  on  a  number 
of  test  cases.  However,  the  “standard”  higher-order 
NURBS  elements  of  type  Qk,  /c  >  2,  produced  poor 
displacements  and  especially  poor  stresses,  whereas 
the  B  and  F  elements  produced  accurate  displace¬ 
ments  and  stresses  in  all  cases.  An  initial  applica¬ 
tion  to  a  problem  of  plasticity  was  presented  but 
the  subject  was  not  systematically  pursued. 

The  present  paper  picks  up  where  Elguedj 
et  al.  [1]  left  off  and  investigates  the  behavior 
of  higher-order  NURBS  elements  on  problems  of 
large-deformation  elastoplasticity.  We  utilize  a 
hyperelastic-based  elastoplastic  theory,  presented 
in  Simo  and  Hughes  [4],  that  has  been  employed 
in  many  investigations  of  numerical  methods.  The 
theory  is  based  on  the  now  classical  F  = 
multiplicative  split  of  the  deformation  gradient  into 
elastic  and  plastic  parts.  The  deformation  gradient 
used  in  the  numerics  is  the  F  mentioned  previously. 
Based  on  numerical  studies  of  benchmark  problems, 
our  experiences  in  large-deformation  elastoplastic¬ 
ity  are  somewhat  different  than  those  for  nonlinear 
elasticity.  Displacements  for  standard  higher-order 
NURBS  elements  are  sometimes  as  good  as  for  F 
elements.  However,  stresses  are  a  much  more  dis¬ 
cerning  measure.  Stresses  are  very  good  for  the  F 
elements  but,  typically,  very  bad  for  the  standard 
elements.  Stresses  invariably  are  the  critical  differ¬ 
entiating  measure  in  problems  of  large- deformat  ion 
plasticity. 

An  outline  of  the  paper  follows:  in  Section  2 
we  very  briefly  introduce  some  ideas  from  NURBS- 
based  IGA.  For  more  extensive  treatments  we  rec¬ 
ommend  Hughes  et  al.  [5],  Cottrell  et  al.  [6]  and, 
in  particular,  Elguedj  et  al.  [1],  in  which  many  re¬ 
sults  pertinent  to  this  work  are  presented.  The 
F  =  F^F^  formulation  of  plasticity  is  briefly  re¬ 
viewed  in  Section  3  and  the  strain  energy  function, 
Mises-Huber  yield  surface,  nonlinear  isotropic  hard¬ 
ening  law,  and  associative  flow  rule  are  defined.  The 
F  method  is  also  described  in  Section  3,  along  with 
the  constitutive  algorithm  and  consistent  elasto¬ 
plastic  moduli.  Numerical  examples  are  presented 
in  Section  4.  The  evaluations  of  the  procedures 
and  comparisons  with  the  Ql/PO  element  results  of 
Simo  and  Armero  [7],  the  mixed  method  elements 
of  Taylor  [8],  and  the  standard  NURBS  elements 
(be.,  those  not  using  the  F  methodology),  are  the 


main  contributions  of  this  paper.  Several  themes 
emerge: 

1.  As  indicated  previously,  displacements  can  be 
very  misleading  for  higher-order  elements.  In 
particular,  the  idea  that  larger  displacements 
are  better  does  not  hold  up  under  scrutiny. 
Standard  higher-order  NURBS  elements  of 
type  Qk  do  not  lock  and  often  produce  dis¬ 
placements  as  large  as  those  for  corresponding 
F  elements,  but  their  stress  distributions  are 
invariably  oscillatory,  revealing  serious  defects. 

2.  Plots  of  reaction  forces  versus  displacements 
can  also  be  misleading. 

3.  Contour  plots  of  stresses  immediately  reveal 
problems  with  element  formulations,  as  will  be 
evident  from  the  results  presented  in  the  se¬ 
quel. 

4.  It  is  important  to  study  refined  finite  ele¬ 
ment  meshes  because  ill-conceived  elements 
may  produce  results  on  refined  meshes  that  are 
no  better  than  for  coarse  meshes.  Obviously, 
the  opposite  holds  true  for  well-conceived  ele¬ 
ments.  We  note  that  most  studies  presented  in 
the  literature  heretofore  only  present  solutions 
on  very  coarse  meshes. 

Our  conclusions  and  our  suggestions  for  further 
research  are  presented  in  Section  5. 

2.  NURBS-based  isogeometric  analysis 
2.1.  Basics 

We  briefly  review  the  concept  of  NURBS-based 
IGA  first  presented  in  Hughes  et  al.  [5],  Cottrell 
et  al.  [6],  where  a  detailed  account  may  be  found. 
NURBS  are  a  generalization  of  B-splines  and  stan¬ 
dard  in  CAD  and  computer  graphics  for  geometry 
modeling  (see,  e.g.^  Cohen  et  al.  [9],  Piegl  and  Tiller 
[10],  Farin  [11]). 

B-splines  are  piecewise  polynomial  functions  with 
a  prescribed  degree  of  continuity.  Univariate  B- 
spline  basis  functions  are  constructed  from  a  knot 
vector,  a  set  of  coordinates  in  parametric  space,  E  = 
{6,6^  •  •  •  :Cn+p+i},  where  G  M  is  the  knot, 
i  is  the  knot  index,  i  =  1,2,  ...,n  +  p+  l,  p  is 
the  polynomial  order,  and  n  is  the  number  of  basis 
functions.  More  than  one  knot  can  be  placed  at  the 
same  location  in  the  parametric  space.  If  m  is  the 
multiplicity  of  a  given  knot,  the  functions  are 
continuous  at  that  location.  If  the  knots  are  equally 
spaced,  the  knot  vector  is  said  to  be  uniform.  A 
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knot  vector  is  referred  to  as  open  if  its  first  and  last 
knots  have  multiplicity  p  +  1.  This  results  in  the 
basis  being  interpolatory  at  the  endpoints  of  the 
interval. 

B-spline  basis  functions  for  a  given  order  p,  are 
defined  recursively  in  the  parametric  space  by  way 
of  the  knot  vector  S.  Beginning  with  piecewise  con¬ 
stants  (p  =  0)  we  have 


Ni, oiO 


1 

0  otherwise. 


For  p  =  1,2,3,...,  the  basis  is  defined  by  the  Cox- 
de  Boor  recursion  formula  (see,  e.g.,  Cohen  et  al. 

[9]): 

NiA^)  = 

Si+p  si 

+  ^i+hP-iiO-  (2) 

Si+p+1  si+1 

Let  ds  denote  the  number  of  spatial  dimensions. 
A  B-spline  curve  in  is  defined  as  follows: 

n 

c(e)  =  ^PiV,p(e,  (3) 

i=l 


where  G  denotes  control  point  i. 

The  univariate  B-spline  concept  can  be  extended 
to  multiple  dimensions  with  the  use  of  tensor  prod¬ 
ucts,  but  representation  of  many  desired  shapes 
of  engineering  interests,  such  as  conic  sections,  re¬ 
quire  further  generalization.  Non-Uniform  Rational 
B-Splines  (NURBS),  rational  projections  of  higher 
dimensional  B-splines,  can  be  introduced  for  this 
purpose  and  consequently  share  many  of  the  same 
properties  as  B-splines.  More  details  can  be  found 
in  Cohen  et  al.  [9],  Rogers  [12]  and  Piegl  and  Tiller 
[10],  as  well  as  Cottrell  et  al.  [6]. 


knot  multiplicity.  The  main  difference  in  IGA  when 
comparing  refinement  strategies  with  finite  element 
ones  is  that  /i-and  p-refinement  do  not  commute. 
The  flexibility  of  knot  insertion  and  order  eleva¬ 
tion  allows  us  to  introduce  /^-refinement  in  which 
order  and  continuity  of  the  basis  functions  can  be 
simultaneously  increased.  This  can  be  attained  by 
performing  order  elevation  on  the  coarse  geometric 
mesh  followed  by  knot  insertion  up  to  the  desired 
mesh  refinement.  It  is  important  to  note  that  k- 
refinement  will  produce  maximum  continuity  on  a 
patch  if  the  coarsest  mesh  is  comprised  of  a  single 
element.  If  the  initial  mesh  comprises  constraints 
on  continuity  across  element  boundaries,  these  will 
exist  on  all  meshes.  For  more  insight  on  mesh  gen¬ 
eration  and  refinement  in  IGA  see  Gottrell  et  al.  [6, 
Chapter  2]. 


3.  Nearly  incompressible  nonlinear  elasticity 
and  plasticity 

3.1.  Elastoplasticity  at  finite  strains 

In  this  section,  we  briefly  review  the  basics  of 
elastoplastic  constitutive  models  at  finite  strains 
based  on  the  notion  of  an  intermediate  stress  free 
configuration.  More  details  can  be  found  for  ex¬ 
ample  in  Simo  [13,  14]  and  Simo  and  Hughes  [4, 
Chapter  9]. 

3.1.1.  Kinematic  preliminaries 

Let  S  C  be  the  reference  configuration  of 
a  body  and  B'  C  be  its  deformed  configura¬ 
tion.  We  introduce  a  mapping  0  :  S  S',  the 
deformation,  that  takes  a  point  X  G  S  to  a  point 
X  =  0(X)  G  S'.  The  deformation  gradient  is  the 
second-order  tensor  defined  by 


2.2.  k-refinement 

Analogues  of  finite  element  /i-and  p-refinement 
are  available  in  IGA.  /^-refinement  is  termed  knot 
insertion  and  consists  in  adding  new  knots  in  the 
knot  vectors.  Knot  insertion  does  not  change  the 
geometric  modeling  and  preserves  continuity  as 
long  as  new  knots  are  not  already  present  in  the 
knot  vector,  p-refinement  is  termed  order  eleva¬ 
tion  and  consists  in  increasing  the  polynomial  or¬ 
der  of  the  basis  functions.  As  with  knot  inser¬ 
tion,  neither  the  geometry  nor  the  parametrization 
is  changed  during  the  process.  Moreover,  continu¬ 
ity  is  preserved  at  element  boundaries  by  increasing 


F  =  V^0(X)  =  ^|^.  (4) 

The  considered  formulation  is  based  on  two  hy¬ 
potheses: 

1.  The  introduction  of  an  intermediate  local  con¬ 
figuration,  relative  to  which  the  elastic  re¬ 
sponse  is  characterized,  which  leads  to  a  mul¬ 
tiplicative  decomposition  of  the  deformation 
gradient: 

F  =  F^FP.  (5) 

2.  The  principle  of  maximum  plastic  dissipation. 
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Following  Eq.  (5),  we  can  introduce  elastic  Eule- 
rian  tensors  and  plastic  Lagrangian  ones  as: 

=  and  e®  =  ^(1  -  (b®)-^),  (6) 

EP  =  ^{CP  -1).  (7) 

Using  the  previous  equations,  we  can  write: 

be  ^  FC^-^F^.  (8) 

Using  the  definition  of  the  total  rate  of  deformation 
tensor,  we  can  express  the  total  rate  of  the  elastic 
left  Cauchy  Green  tensor  as  a  Lie  derivative  by  the 
following  equation: 

=  FCP“^F^.  (9) 

Within  the  context  of  the  infinitesimal  theory, 
the  small  strain  tensor  s  is  decomposed  into  volu¬ 
metric  and  deviatoric  parts  using  an  additive  split. 
In  the  nonlinear  theory,  the  decomposition  becomes 
multiplicative: 

F  =  (10) 

where 

detF  =  J  =  detF'^“  and  detF'’®''  =  1,  (11) 

which  leads  to: 

Fdev  ^  J-l/3p  pdil  ^  jl/3j_  (^2) 

3.1.2.  J2  flow  theory  at  finite  strains 
We  consider  a  J2  how  theory  with  isotropic  hard¬ 
ening.  We  assume  that  the  stress  response  is 
isotropic  and  that  the  plastic  flow  is  isochoric,  that 
is: 

detF^  =  detC^  =  1  and 
J  =  detF  =  detFU  (13) 

Using  the  assumption  of  isotropy  and  the  in¬ 
termediate  stress  free  configuration,  the  stress  re¬ 
sponse  can  be  characterized  using  a  nonlinear  elas¬ 
tic  strain  energy  function: 

w  =  u{r)^w{h^),  with 
b®  =  (J®)-2/3b®. 


In  the  remaining,  we  will  use  the  following  expres¬ 
sions  for  the  isochoric  and  volumetric  parts  of  the 


strain  energy  function: 

U(J)  =  1.  (1(. 

r2-l)-lnJ^, 

(15) 

W{h)  = 

(tr[b]  -  3)  , 

(16) 

where  n  is  interpreted  as  the  bulk  modulus  and  /i 
the  shear  modulus  in  infinitesimal  deformations. 

We  consider  the  classical  Mises-Huber  yield  con¬ 
dition,  formulated  in  terms  of  the  Kirchhoff  stress 
tensor  as: 

f{T,a)  =  ||dev[T]||  -  y|fc(a)  <  0.  (17) 

We  will  consider  in  this  paper  a  nonlinear  isotropic 
hardening  law  using  the  following  expression 
(see,  e.g.^  Simo  and  Hughes  [4]): 

k{a)  =  (Jo  +  (cToo  -  cro)[l  -  exp(-^(a)] 

+  Ka^  with  ^  >  0,  (18) 

where  ao  is  the  initial  flow  stress,  (Too  is  the  sat¬ 
uration  flow  stress,  5  is  the  saturation  exponent, 
K  is  the  linear  hardening  coefficient,  and  a  is  the 
equivalent  plastic  strain. 

As  in  the  infinitesimal  theory,  the  associative  flow 
rule  is  uniquely  determined  by  the  principle  of  max¬ 
imum  plastic  dissipation  given  the  strain  energy 
function  and  the  yield  condition: 

i„b®  =  -?7tr[b®]n,  (19) 

n  =  s/||s||,  (20) 

s  =  dev[r].  (21) 

The  evolution  of  the  hardening  variable  is  the 
same  as  in  the  infinitesimal  theory: 


where  7  is  the  consistency  parameter  subject  to  the 
Kuhn- Tucker  conditions: 

7  >  0,  /(r,  a)  <  0,  7/(t,  a)  =  0,  (23) 


(14) 
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with  the  consistency  condition 
7/(t,q;)  =  0. 


(24) 


3.1.3.  Integration  algorithm 

The  integration  algorithm  of  the  previous  formu¬ 
lation  can  be  obtained  using  the  general  procedure 
presented  in  Simo  and  Hughes  [4].  We  proceed  as 
follows: 

•  The  equations  are  pull-backed  to  the  material 
description  using  the  deformation  gradient. 

•  The  time  stepping  algorithm  is  performed  in 
the  material  description.  In  our  case  we  use 
backward  Euler. 

•  The  discrete  evolution  equations  are  pushed 
forward  to  the  spatial  description  using  the  de¬ 
formation  gradient. 

Using  the  obtained  equations,  the  return  mapping 
algorithm  can  be  constructed,  the  problem  being 
strain  driven.  This  algorithm  is  summarized  in 
Box  1.  In  addition,  we  present  in  the  appendix 
the  expression  for  the  consistent  elastoplastic  mod¬ 
uli  for  the  proposed  radial  return  algorithm  (see 
Simo  [13,  14])  as  well  as  the  iterative  algorithm  used 
to  solve  the  consistency  condition  for  the  proposed 
nonlinear  isotropic  hardening  law. 

3.2.  F  formulation  for  nearly  ineompressible  non¬ 
linear  elastieity  and  plastieity 

The  F  method  is  a  generalization  of  the  B 
method  (see,  e.g.,  Hughes  [2,  3]),  a  strain-projection 
technique,  to  the  hyperelastic  finite-deformation 
regime.  The  central  idea  is  to  multiplicatively  split 
the  deformation  gradient  tensor  F  into  its  devia- 
toric  (volume  preserving)  and  dilatational  (volume 
changing)  parts  as  given  in  Eq.  (10).  This  decompo¬ 
sition  has  been  exploited  previously  by  Flory  [15], 
Hughes  et  al.  [16],  Simo  et  al.  [17],  Simo  and  Tay¬ 
lor  [18]  (within  a  three  field  Hu-Washizu  principle), 
and,  more  recently,  de  Souza  Neto  et  al.  [19,  20]  in 
an  alternative  F  approach. 

We  define  a  modified  deformation  gradient  F  in 
terms  of  F^®^  and  a  modified  dilatational  part  of 
the  deformation  gradient  F^^^: 

F  =  F^i^F^®^,  (25) 

where 

pdii  _  _  P(  (26) 

with  P  a  linear  projection  operator  with  respect  to 
the  I/^-norm  onto  a  space  of  lower-order  and  lower- 
continuity  functions.  This  allows  us  to  obtain  a 


robust  formulation  with  only  displacements  as  pri¬ 
mary  unknowns  in  order  to  solve  large  strain  nearly 
incompressible  problems.  For  piecewise  polyno¬ 
mial  basis  functions  such  as  NURBS,  the  projec¬ 
tion  space  is  taken  as  the  space  one  order  lower 
than  the  displacement  space.  This  introduces  the 
family  of  elements  Qk/Qk-1  where  k  is  the  polyno¬ 
mial  order  of  the  displacement  space.  The  space 
Qk  is  continuous  and  the  space  Qk-1  is 

continuous.  For  further  details  the  interested 
reader  should  refer  to  Elguedj  et  al.  [1]. 

We  employ  here  the  proposed  F  method  in  con¬ 
junction  with  the  previously  presented  elastoplas¬ 
tic  model.  In  practice  this  comes  down  to  using 
the  modified  deformation  gradient  F  instead  of  the 
original  deformation  gradient  F  as  an  input  in  all 
the  elastoplastic  routines.  The  internal  force  vec¬ 
tor  and  consistent  tangent  stiffness  matrix  for  hy¬ 
perelasticity  are  defined  in  Elguedj  et  al.  [1].  For 
the  elastoplastic  theory  used  here,  we  construct  the 
consistent  tangent  stiffness  matrix  by  replacing  the 
tensor  of  elastic  moduli  in  the  formula  derived  in 
Elguedj  et  al.  [1]  by  the  elastic-plastic  moduli  pre¬ 
sented  in  Box  2. 

4.  Examples 

In  all  the  examples  presented  in  the  following 
subsections,  the  material  model  consists  of  a  Neo- 
Hookean  one  for  the  elastic  part,  and  an  associative 
flow  rule  based  on  a  Von  Mises  yield  criterion  with 
nonlinear  isotropic  hardening  following  a  saturation 
law  for  the  plastic  part.  The  nonlinear  isotropic 
hardening  rule  is  defined  from  Eqs.  (17)  and  (18). 
The  strain  energy  function  of  the  Neo-Hookean 
model  is  given  by  Eqs.  (15)  and  (16).  The  same 
material  parameters  are  used  in  all  the  examples 
and  are  given  in  Table  1. 


Shear  modulus  /r 

80.1938  GPa 

Bulk  modulus  n 

164.21  GPa 

Initial  flow  stress  ctq 

450  MPa 

Saturation  flow  stress  (Too 

715  MPa 

Saturation  exponent  5 

16.93 

Linear  hardening  coefficient  K 

129.24  MPa 

Table  1:  Elastic-plastic  material  properties. 


4.1.  Cook’s  membrane 

We  consider  as  a  first  example  the  elastic-plastic 
extension  of  the  well  known  Cook’s  membrane. 
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1.  Update  current  configuration 


O  0n  —  ^n(^)  +  ^^[^^(X)], 

f  .  —  1  _j_ 
tn+i  -  i  H- 

Fn+l  —  fn+lFn- 

2.  Compute  elastic  predictor 

fn+i  =  (det[f„+i])“^^^f„+i, 

Ce  trial  _  f  Ce  f  T 

^n+1  , 

s*„”f  =  Aidev[b-^r']. 

3.  Check  for  plastic  loading 

Compute  II  -  ^K{a^). 

if  fn+i  —  0  elastic  state  =>  save  and  exit, 
else  perform  return  mapping. 

4.  Return  mapping 

p  =  , 

Compute  A7  such  that  /(A7)  =  ||s^”“*||  —  f^K{an  +  IA7)  —  2jlA'y  =  0, 

rt  _  f^trial  /\\^trial\\ 

“  —  ®n+l  /  IPn+1  II  5 

Sn+i  =  -  2/2A7n  , 


^n+1  — 

5.  Stress  update 


Jn-\-l  —  det[F77_|_l]  =>  Pn+l  —  U  (t/n+l)  —  2  ('^n+1  ^)/ '^n+l? 

^n+1  ~  'A^+lPn+ll  “1“  S77_|_i. 

6.  Update  intermediate  configuration 

hl+i  =  (1/m)s„+i  +  (/2/m)1. 


Box  1:  Return-mapping  algorithm  for  J2— flow  theory  with  nonlinear  isotropic  hardening. 


This  example  has  been  studied  by  many  researchers 
to  validate  new  methods  to  treat  incompressibility 
under  combined  bending  and  shear  in  both  small 
and  large  deformation  cases.  The  extension  to  large 
strain  plasticity  was  proposed  by  Simo  and  Armero 
[7]  and  later  studied  by  Armero  and  Glaser  [21] 
and  Ramesh  and  Maniatty  [22].  A  tapered  panel 
is  clamped  on  one  side  and  subjected  to  a  uniform 
vertical  shear  load  on  the  opposite  side.  The  other 
sides  are  traction  free.  The  geometry,  loading  and 
boundary  conditions  are  given  in  Figure  1.  The 
load  value  is  chosen  to  be  F  =  bN/mm  and  is  im¬ 
posed  in  ten  equal  steps.  Plain  strain  conditions 
are  assumed  to  apply. 

The  results,  shown  in  Figure  2  for  the  vertical  dis¬ 


placement  of  the  top  right  point  versus  number  of 
elements  per  side,  compare  well  to  those  obtained 
in  Simo  and  Armero  [7]  and  Armero  and  Glaser 
[21].  Similar  convergence  behavior,  when  increas¬ 
ing  the  order  of  approximation  with  the  proposed  F 
method,  to  those  obtained  in  the  nearly  incompress¬ 
ible  linear  and  nonlinear  elastic  cases  (see  Elguedj 
et  al.  [1])  is  obtained.  In  Figure  3,  we  show  the 
ayy  component  of  the  Cauchy  stress  tensor  for  the 
4x4  element  mesh  with  linear,  quadratic  and  cu¬ 
bic  NURBS  with  and  without  F  .  We  can  see  stress 
oscillation  patterns  for  the  non-F  cases,  even  for 
the  cubic  NURBS  that  showed  reasonable  conver¬ 
gence  for  the  top  right  corner  displacement.  This 
behavior  is  also  similar  to  what  was  observed  by 
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Figure  2:  Elastic-plastic  Cook’s  membrane.  Vertical  displacement  of  top  right  corner  versus  number  of 
elements  per  edge  with  and  without  F  for  various  NURBS  orders  obtained  from  /c-refinement. 


Figure  1:  Geometry,  loading,  boundary  conditions, 
material  parameters  and  quantity  of  interest  for  the 
plane  strain  Cook’s  membrane. 


the  authors  in  the  nonlinear  elastic  incompressible 
case  in  [1].  Clear  locking  behavior  may  be  noted 
for  Q1  in  Figure  2,  but  not  Q2  and  Q3.  However, 
stress  oscillations  which  are  especially  apparent  for 
Q2,  are  initial  indications  of  deficiencies  of  non-F 
elements. 

4.2.  Plane  strain  necking  of  a  bar 

In  the  second  example,  we  consider  a  plane-strain 
bar  subjected  to  uniform  extension.  This  is  a  stan¬ 
dard  test  problem  that  was  also  analyzed  by  Simo 
and  Armero  [7],  Armero  and  Glaser  [21],  Dvorkin 
and  Assanelli  [23],  de  Souza  Neto  et  al.  [19,  20], 
agelet  de  Saracibar  et  al.  [24]  and  Taylor  [8].  The 
bar  has  a  length  of  53.334mm  and  a  width  of 
12.826mm,  symmetry  conditions  are  used  and  only 
one  fourth  of  the  domain  is  modeled  in  the  numeri¬ 
cal  simulation  as  presented  in  Figure  4.  In  order  to 
control  the  location  of  the  necking,  the  dimension 
at  the  center  of  the  bar  is  reduced  to  0.982  of  the 
top  width.  A  total  displacement  of  5mm  is  applied 
on  the  top  surface.  Numerical  simulations  are  per¬ 
formed  with  meshes  of  50  and  200  elements  with 
refinement  near  the  necking  area  (see  Figure  4). 

In  Figures  5  and  6  we  plot  the  ayy  component  of 
the  Cauchy  stress  tensor  and  the  equivalent  plastic 
strain  on  the  deformed  configuration  at  maximum 
load  for  linear,  quadratic  and  cubic  cases  with  and 
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Figure  3:  Elastic-plastic  Cook’s  membrane,  coarse 
component  at  maximum  load  for  (a)  Ql,  (c)  Q2,  (e) 


Sigma„vv  / 
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;h  with  four  elements  per  side,  ayy  Cauchy  stress 
,  (b),  F  Ql/QO  (d)  F  Q2/Q1  (f)  F  Q3/Q2. 
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imposed  displacement 
of  5mm 


due  to  symmetry,  only 
1/4  is  modelled 


Uy  =0 


200  elements 


Figure  4:  Plane  strain  necking:  problem  setup  and  quarter  meshes  used  for  the  computation 


without  F  for  the  50  element  case.  All  results  are 
shown  with  different  scales  adapted  to  the  minimum 
and  maximum  values  obtained  in  each  case.  For  the 
non-F  cases,  we  can  clearly  see  oscillation  patterns 
due  to  locking  for  the  stress.  Such  patterns  can¬ 
not  be  seen  in  the  F  cases  even  for  the  linear  case 
with  this  relatively  coarse  mesh.  Although  we  ob¬ 
serve  stress  oscillations,  the  deformed  shapes  of  the 
specimen  look  very  similar  between  the  F  and  non- 
F  cases  for  quadratic  and  cubic  functions.  From 
Figure  5  we  observe  that  the  minimum  and  max¬ 
imum  values  of  the  ayy  component  of  the  Cauchy 
stress  tensor  for  the  piecewise  linear  case  are  every 
different  from  the  values  for  the  quadratic  and  cu¬ 
bic  cases,  for  both  the  F  and  non-F  cases.  These 
values  are  quite  close  to  each  other  when  quadratic 
and  cubic  F  cases  are  compared.  Similar  observa¬ 
tions  can  be  made  for  the  equivalent  plastic  strain 
in  Figure  6. 

In  Figures  7  and  8  we  plot  the  horizontal  displace¬ 
ment  at  the  bottom  right  corner  of  the  specimen 
(z.e.,  where  necking  occurs)  and  the  reaction  force 
versus  the  imposed  vertical  displacement.  Results 
are  plotted  with  dashed  lines  for  non-F  cases  and 
solid  lines  for  F  cases.  We  can  see  that  F  Ql/QO 
provides  reasonable  results  but  that  the  mesh  is  too 
coarse  to  achieve  the  results  of  Q2/Q1  and  Q3/Q2. 


As  for  the  Cook’s  membrane  we  can  note  that  for 
the  non-F  Q2  and  Q3  cases,  although  stress  oscil¬ 
lations  are  observed,  global  parameters  such  as  dis¬ 
placement  and  reaction  force  are  close  to  the  correct 
values. 

In  Figures  9  to  12  we  plot  similar  results  as  pre¬ 
viously  but  for  the  200  element  case.  All  results  are 
shown  with  different  scales  adapted  to  the  minimum 
and  maximum  values  obtained  in  each  case.  We  can 
still  observe  oscillations  for  the  stress  in  the  non-F 
cases.  As  with  the  coarse  mesh,  the  minimum  and 
maximum  values  for  both  quantities  for  the  Q1  case 
are  very  different  from  the  other  cases.  However,  we 
can  see  that  for  Ql/ QO  the  maximum  and  minimum 
values  are  converging  to  the  ones  obtained  in  the 
Q2/Q1  and  Q3/Q2  cases.  For  the  necking  displace¬ 
ment  and  the  reaction  force,  we  also  plotted  results 
obtained  using  mixed  finite  elements  Ql/PO  and 
Q2/P1  quadrilaterals  taken  from  Taylor  [8].  We 
can  see  that  the  results  for  Ql/PO  are  identical  to 
the  ones  obtained  with  Ql/QO  and  that  results  for 
Q2/P1  are  very  close  to  Q2/Q1.  We  note  that  the 
results  for  F  Q3/Q2  and  Q2/Q1  are  indistinguish¬ 
able  (results  are  identical  up  to  the  fifth  digit  for 
both  quantities).  We  also  wish  to  emphasize  that, 
in  our  opinion,  the  necking  displacement  and  reac¬ 
tion  results  may  be  misleading  in  that  even  when 
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Figure  5:  Plane  strain  necking,  coarse  mesh  with  50  elements,  (jyy  Cauchy  stress  component  for  a  vertical 
imposed  displacement  of  5mm  for:  (a)  Ql,  (b)  Q2,  (c)  Q3,  (d)  F  Ql/QO,  (e)  F  Q2/Q1,  (f)  F  Q3/Q2. 
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Figure  6;  Plane  strain  necking,  coarse  mesh  with  50  elements.  Equivalent  plastic  strain  a  for  a  vertical 
imposed  displacement  of  5mm  for:  (a)  Ql,  (b)  Q2,  (c)  Q3,  (d)  F  Ql/QO,  (e)  F  Q2/Q1,  (f)  F  Q3/Q2. 
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Figure  7:  Plane  strain  necking,  coarse  mesh  with  50  elements.  Necking  displacement  versus  imposed  dis¬ 
placement. 


imposed  vertical  displacement  /  mm 


Figure  8:  Plane  strain  necking,  coarse  mesh  with  50  elements.  Half  reaction  force  versus  imposed  displace¬ 
ment. 
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stresses  are  clearly  inaccurate,  such  as  for  the  cases 
Q2  and  Q3,  the  necking  and  reaction  results  are 
very  similar  to  those  for  F  Q2/Q1  and  F  Q3/Q2. 

4.3.  Three  dimensional  neeking 

The  final  example  is  the  three-dimensional  ex¬ 
tension  of  the  previous  one.  It  has  been  studied  by 
Simo  and  Armero  [7],  de  Souza  Neto  et  ah  [19,  20] 
and  agelet  de  Saracibar  et  ah  [24].  A  cylindrical 
bar  with  the  same  dimensions  as  in  the  2D  case  is 
considered  and,  due  to  symmetry,  only  one  eighth 
of  the  bar  is  modeled  in  the  numerical  simulation, 
see  Figure  13.  In  order  to  trigger  the  necking,  we 
also  reduce  the  radius  of  the  bar  in  the  center  of  the 
specimen  to  0.982  of  the  top  radius.  We  used  two 
meshes  of  4  x  4  x  10  and  8  x  8  x  20  elements  with  re¬ 
finement  near  the  necking  area  (see  Figure  13),  and 
the  loading  is  imposed  as  a  vertical  displacement  on 
the  top  surface  up  to  9mm.  Note  here  that  other 
authors  imposed  a  displacement  of  7mm  only.  The 
maximum  imposed  displacement  was  increased  here 
in  order  to  fully  distinguish  results  between  F  and 
non-F  cases  as  will  be  seen  subsequently. 

We  investigate  for  this  example  the  behavior  of 
the  linearized  operator  described  previously.  We 
see  in  Table  2  the  evolution  of  the  relative  Euclid¬ 
ian  norm  of  the  residual  over  the  iterations  for  the 
Newton  step  corresponding  to  an  imposed  displace¬ 
ment  of  7mm  for  the  coarse  mesh  with  Q2,  Q3,  F 
Q2/Q1  and  F  Q3/Q2.  We  can  observe  that  conver¬ 
gence  is  slightly  better  for  F  cases  when  compared 
with  the  non-F  cases. 

In  Figures  14  and  15,  we  plot  for  the  coarse  mesh 
the  azz  component  of  the  Cauchy  stress  tensor  and 
the  equivalent  plastic  strain  on  the  deformed  con¬ 
figuration  for  an  imposed  displacement  of  7mm  for 
quadratic  and  cubic  NURBS  with  and  without  F. 
Although  results  were  computed  on  one  eighth  of 
the  specimen,  the  results  were  reflected  for  bet¬ 
ter  understanding.  The  scales  are  adapted  for  each 
case  with  different  minimum  and  maximum  values. 
Again  we  observe  oscillation  patterns  for  the  non-F 
cases  and  smooth  fields  for  the  F  cases  when  look¬ 
ing  at  the  azz  component  of  the  Cauchy  stress  ten¬ 
sor.  Such  oscillation  patterns  are  not  observed  for 
the  equivalent  plastic  strain  but  the  minimum  and 
maximum  values  obtained  for  Q2  are  quite  different 
from  the  three  other  cases. 

In  Figures  16  and  17  we  plot  the  necking  dis¬ 
placement  and  reaction  force  versus  applied  load 
for  the  coarse  mesh.  The  results  are  shown  for 
quadratic  and  cubic  NURBS  with  and  without  F 


and  for  the  reaction  force  we  also  plotted  results 
obtained  with  the  mixed  finite  element  Ql/PO  ex¬ 
tracted  from  Simo  and  Armero  [7]  for  comparison. 
We  can  see  that  up  to  6mm  loading,  the  Q2  case 
produces  a  reasonable  result.  This  is  also  the  case 
for  Q3  up  to  7mm,  which  is  the  reason  that  mo¬ 
tivated  the  increase  of  the  maximum  imposed  dis¬ 
placement.  When  the  load  continues  to  increase, 
the  solution  becomes  stiffer  for  the  non-F  cases  and 
deviates  from  the  reference  one.  We  can  see  that 
both  F  cases  are  close  to  each  other.  A  similar 
picture  is  observed  when  looking  at  the  reaction 
force.  However,  we  observe  that  Q3  is  very  close  to 
Q2/Q1  up  to  the  maximum  load  and  that  the  dif¬ 
ference  between  Q2/Q1  and  Q3/Q2  is  a  little  more 
pronounced  than  for  the  necking  displacement.  We 
can  see  that  up  to  6mm  both  F  cases  compare  well 
with  Ql/PO  that  was  obtained  with  a  much  finer 
mesh. 

We  plot  similar  results  obtained  with  quadratic 
NURBS  on  the  fine  mesh  in  Figures  18  to  20.  For 
the  azz  component  of  the  Cauchy  stress  tensor,  we 
still  observe  oscillations  for  Q2  but  not  for  Q2/Q1. 
Interestingly,  although  necking  displacement  and 
reaction  force  results  are  quite  good  for  Q2  at  6mm 
with  the  fine  mesh  and  close  to  Q3  results  obtained 
on  the  coarse  mesh,  the  stress  oscillations  produce 
minimum  and  maximum  stress  values  that  are  very 
similar  to  Q2  on  the  coarse  mesh.  This  indicates 
that  although  the  mesh  is  finer,  stress  oscillations 
are  not  reduced  with  quadratic  NURBS.  Results 
obtained  for  the  equivalent  plastic  strain  are  very 
similar  to  the  ones  obtained  for  the  coarse  mesh  and 
Q3/Q2,  although  the  maximum  value  for  Q2/Q1 
on  the  fine  mesh  is  10%  higher.  When  looking 
at  necking  displacement  and  reaction  force  curves 
for  Q2/Q1,  we  can  see  that  up  to  7mm  the  results 
are  almost  identical  to  Q3/Q2  on  the  coarse  mesh. 
When  the  load  increases  further,  the  response  is  a 
little  softer  with  the  finer  mesh,  but  the  first  mesh  is 
probably  too  coarse  in  the  necking  area  to  capture 
the  complete  behavior  at  such  a  high  loading. 

Finally,  we  plot  in  Figures  21  and  22  the  de¬ 
formed  meshes  for  an  imposed  displacement  of  7mm 
and  9mm  for  F  and  non-F  cubic  NURBS  with  the 
coarse  mesh  and  quadratic  NURBS  with  the  fine 
mesh.  We  can  see  for  an  imposed  displacement  of 
7mm,  as  was  observed  on  the  necking  displacement 
curves,  that  all  cases  shown  are  very  similar.  Al¬ 
though  stress  oscillations  are  present  for  the  non 
F  cases,  the  deformed  configuration  and  meshes  do 
not  exhibit  unusual  behavior.  For  an  imposed  dis- 
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Figure  9:  Plane  strain  necking,  fine  mesh  with  200  elements,  (jyy  Cauchy  stress  component  for  a  vertical 
imposed  displacement  of  5mm  for:  (a)  Ql,  (b)  Q2,  (c)  Q3,  (d)  F  Ql/QO,  (e)  F  Q2/Q1,  (f)  F  Q3/Q2. 
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Figure  10;  Plane  strain  necking,  fine  mesh  with  200  elements.  Equivalent  plastic  strain  a  for  a  vertical 
imposed  displacement  of  5mm  for:  (a)  Ql,  (b)  Q2,  (c)  Q3,  (d)  F  Ql/QO,  (e)  F  Q2/Q1,  (f)  F  Q3/Q2. 


15 
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Figure  11:  Plane  strain  necking,  fine  mesh  with  200  elements.  Necking  displacement  versus  imposed  dis¬ 
placement.  Mixed  finite  elements  are  taken  from  Taylor  [8]. 


Figure  12:  Plane  strain  necking,  fine  mesh  with  200  elements.  Half  reaction  force  versus  imposed  displace¬ 
ment.  Mixed  finite  elements  are  taken  from  Taylor  [8]. 
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Iteration  number 

Relative  norm 

of  the  residual 

Q2 

Q3 

F  Q2/Q1 

F  Q3/Q2 

1 

1.000000x10° 

1.000000x10° 

1.000000x10° 

1.000000x10° 

2 

0.062626x10° 

0.0550812x10° 

0.064762x10° 

0.052881x10° 

3 

0.952021x10-3 

0.0047900x10° 

0.777876x10-3 

0.562725x10-3 

4 

0.621975x10-4 

0.335220x10-3 

0.130120x10-3 

0.955583x10-° 

5 

0.865509x10-'^ 

0.467890x10-4 

0.246640x10-3 

0.183784x10-3 

6 

0.129028x10-° 

0.756419x10-^ 

0.171958x10-4° 

0.229679x10-4° 

7 

0.307675x10-“ 

0.593532x10-4° 

Table  2:  3D  necking.  Evolution  of  the  relative  norm  of  the  residual  during  the  Newton  step  with  imposed 
displacement  of  7mm  for  the  coarse  mesh  with  Q2,  Q3,  F  Q2/Q1  and  F  Q3/Q2. 


placement  of  9mm,  we  can  see  that  the  deformation 
is  too  localized  in  the  first  layer  of  elements  and  that 


17 


it  is  difficult  to  obtain  an  accurate  solution.  How¬ 
ever,  we  can  see,  as  was  observed  on  the  necking 
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Figure  14:  3D  necking,  coarse  mesh,  azz  component  of  the  Cauchy  stress  tensor  on  the  deformed  configu¬ 
ration  for  a  vertical  imposed  displacement  of  7mm.  (a)  Q2,  (b)  Q3,  (c)  F  Q2/Q1,  (d)  F  Q3/Q2. 
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Figure  15:  3D  necking,  coarse  mesh.  Equivalent  plastic  strain  ol  on  the  deformed  configuration  for  a  vertical 
imposed  displacement  of  7mm.  (a)  Q2,  (b)  Q3,  (c)  F  Q2/Q1,  (d)  F  Q3/Q2. 
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imposed  vertical  displacement/mm 


Figure  16:  3D  necking,  coarse  mesh.  Necking  displacement  versus  imposed  displacement. 


Figure  17:  3D  necking,  coarse  mesh.  Reaction  force  versus  imposed  displacement.  Ql/PO  taken  from  Simo 
and  Armero  [7],  for  a  calculation  with  six  times  as  many  elements  as  the  present  meshes. 
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Figure  18:  3D  necking,  fine  mesh,  azz  component  of  the  Cauchy  stress  tensor  and  equivalent  plastic  strain 
on  the  deformed  configuration  for  a  vertical  imposed  displacement  of  7mm.  (a)  azz  Q2,(b)  Q2,  (c)  azz 

Q2/Q1,  (d)  a  F  Q2/Q1. 
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Figure  19:  3D  necking,  fine  mesh.  Necking  displacement  versus  imposed  displacement. 


Figure  20:  3D  necking,  fine  mesh.  Reaction  force  versus  imposed  displacement.  Ql/PO  taken  from  Simo 
and  Armero  [7] 
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displacement  curves,  that  the  non-F  cases  do  pro¬ 
duce  a  stiffer  response  compared  to  the  F  ones. 

5.  Conclusion 

In  this  paper  we  have  presented  an  evalua¬ 
tion  of  the  performance  of  NURBS-based  Isoge¬ 
ometric  Analysis  elements  on  problems  of  large- 
deformation  elastoplasticity.  The  elastoplastic  con¬ 
stitutive  model  utilizes  the  multiplicative  splitting 
of  the  deformation  gradient  into  elastic  and  plastic 
parts,  F  =  F^F^.  The  element  technology  is  based 
on  the  F  method  in  which  the  dilatational  compo¬ 
nent  is  projected  on  a  lower-order  space  than  the 
one  used  for  the  displacements.  The  F  method, 
which  is  derived  from  a  modified  minimum  poten¬ 
tial  energy  principle,  is  amenable  to  exact  and  ex¬ 
plicit  linearization,  and  produces  a  symmetric  con¬ 
sistent  tangent  matrix.  Normally,  the  reason  to  em¬ 
ploy  an  F  approach  is  thought  to  be  to  circumvent 
mesh  locking,  which  is  manifest  for  the  standard  Q1 
element.  However,  mesh  locking  is  not  the  issue  for 
higher-order  elements  in  plasticity  as  our  numerical 
results  revealed.  We  studied  the  standard  NURBS 
elements  (be.,  no  F)  Ql,  Q2,  Q3,  and  Q4,  and  com¬ 
pared  them  with  the  F  elements  Ql/QO,  Q2/Q1, 
Q3/Q2,  and  Q4/Q3,  and  with  the  Ql/PO  element 
results  of  Simo  and  Armero  [7],  and  the  mixed  el¬ 
ements  of  Taylor  [8].  Our  computations  provided 
the  following  insights  into  the  evaluation  of  higher- 
order  NURBS  elements  in  plasticity: 

1.  Point  displacements  are  a  misleading  metric. 

2.  Graphs  of  reaction  forces  versus  displacements 
may  also  be  misleading. 

3.  Stress  distributions,  in  the  form  of  contour 
plots,  expose  the  deficiencies  of  ill-conceived  el¬ 
ements. 

4.  Mesh  refinement  studies  are  important  indica¬ 
tors  because  results  for  ill-conceived  elements 
are  often  no  better  for  refined  meshes  than  for 
coarse  meshes. 

In  future  work,  we  recommend  performing  seri¬ 
ous  mesh  refinement  studies  to  converge  stress  and 
equivalent  plastic  strains,  as  well  as  displacements 
and  reaction  forces.  Of  course,  such  studies  will  re¬ 
quire  significant  computational  resources  and  are 
easier  said  than  done.  Nevertheless,  a  database 
of  results  obtained  would  constitute  an  invaluable 
benchmarking  tool  and  would  eliminate  some  of  the 
ambiguity  in  results  that  have  appeared  in  the  lit¬ 
erature.  The  plane  strain  and  3D  necking  problems 


would  seem  appropriate  for  such  an  exercise.  We 
feel  the  so-called  Cooks  membrane  is  not  ideal  for 
this  purpose  because  the  singularity  at  the  reen¬ 
trant  corner  dominates  and  clouds  evaluation  as 
emphasized  to  us  by  Kjell  Magne  Mathisen. 

We  believe  our  F  formulation  utilizing  the 
smooth  NURBS  spaces  Qk/Qk-1  has  been  demon¬ 
strated  herein  to  be  an  effective  approach  to  solving 
large-deformation  plasticity  problems  and  has  been 
shown  to  possess  several  attributes.  However,  it 
has  some  shortcomings  concerning  the  structure  of 
the  consistent  tangent  stiffness  matrix  and  numer¬ 
ical  quadrature  that  should  be  addressed  in  future 
work.  The  tangent  stiffness  exhibits  dense  coupling 
on  NURBS  patches  in  the  case  of  k  >  2,  due  to 
the  continuity  of  the  F  Jacobian  determinant.  In 
Elguedj  et  al.  [1]  we  described  an  effective  equa¬ 
tion  solution  strategy  for  circumventing  this  cou¬ 
pling,  but  this  strategy  is  not  one  that  is  currently 
supported  in  commercial  finite  element  programs. 
Thus,  it  would  be  desirable  to  develop  an  equally  ef¬ 
fective  generalization  of  the  method  that  eliminated 
the  coupling  ab  initio.  We  have  employed  Gauss 
quadrature  rules  of  order  /c  +  1  in  each  coordinate 
direction  on  Bezier  elements  (see,  e.^.,  Scott  et  al. 
[25],  Borden  et  al.  [26]).  This  is  necessary  to  attain 
full  rank  for  Ql/QO,  but  may  amount  to  overkill 
for  higher-order  elements.  It  is  always  desirable  to 
use  quadrature  rules  with  the  minimum  number  of 
points  in  order  to  save  computational  effort.  The 
challenge  is  to  find  minimum-point  rules  that  main¬ 
tain  stability  and  do  not  otherwise  degrade  perfor¬ 
mance.  Efforts  to  develop  efficient  quadrature  rules 
for  higher-order  NURBS  spaces  were  initiated  in 
Hughes  et  al.  [27].  Another  direction  that  is  appeal¬ 
ing,  due  to  its  utter  simplicity,  is  to  use  lower-order 
“reduced  quadrature”.  In  some  academic  circles 
this  is  considered  anathema,  but  the  fact  remains 
it  is  successfully  employed  in  widely  used  commer¬ 
cial  finite  element  programs  and,  in  our  opinion,  is 
certainly  worth  exploring.  Benson  et  al.  [28]  have 
reported  initial  success  with  this  approach. 
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(a)  (b)  (c)  (d) 


Figure  21:  3D  necking,  deformed  meshes  at  u=7mm.  (a)  Q3  coarse  mesh,  (b)  Q2  fine  mesh,  (c)  Q3/Q2 
coarse  mesh,  (d)  Q2/Q1  fine  mesh. 


(a)  (b)  (c)  (d) 

Figure  22:  3D  necking,  deformed  meshes  at  u=9mm.  (a)  Q3  coarse  mesh,  (b)  Q2  fine  mesh,  (c)  Q3/Q2 
coarse  mesh,  (d)  Q2/Q1  fine  mesh. 
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A.  Consistent  elastoplastic  tangent  moduli 

The  algorithm  given  in  Box  1  is  amenable  to  ex¬ 
act  linearization  leading  to  a  closed  form  expression 
for  the  consistent  algorithmic  tangent  moduli  in  the 
finite  strain  range.  The  linearization  can  be  carried 
out  in  closed  form  because  of  the  hyperelastic  na¬ 
ture  of  the  stress  response.  The  complete  deriva¬ 
tion  can  be  found  in  Simo  [13,  14],  the  equations 
are  summarized  in  Box  2. 

B.  Consistency  condition 

Although  the  counterpart  of  the  consistency  con¬ 
dition  can  be  determined  in  closed  form  for  the  case 
of  linear  hardening,  it  requires  an  iterative  algo¬ 
rithm  for  the  general  non-linear  case  considered  in 
this  paper.  The  consistency  condition  is  given  in 
Eq.  (17).  For  the  nonlinear  hardening  law  we  have 
considered,  its  counterpart  is  written  as: 

/(A7)  =  iis*„Ti'ii  -  +  y|A7) 

-  2/IA7  =  0.  (27) 

This  expression  furnishes  a  nonlinear  equation  for 
Ay  which  is  easily  solved  using  an  iterative  method. 
If  the  derivative  of  k{a)  is  easily  computed  in  closed 
form,  a  Newton  iteration  scheme  can  be  employed 
(see  Simo  and  Hughes  [4],  Simo  [14],  Simo  and  Tay¬ 
lor  [29]).  The  algorithm  is  given  in  Box  3  and  is 
guaranteed  to  converge  at  a  quadratic  rate  pro¬ 
vided  that  —k{a)  is  convex,  which  is  the  case  for 
the  model  given  in  Eq.  (18). 
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